clear; clc; close all;
angstrom = 1e-10;  %[m]
a0 = 5.29177e-11;  %[m]
cmn1 = 0.0299792458e12; %[Hz/cm^-1]
foffset = -0.08;     %[GHz] wavemeter offset(-0.16 GHz) + AOM(+0.08GHz)

fnames = dir('data/v4/2020*Rb2*.mat');
Nfile = length(fnames);
dataAll = cell(1, Nfile);
disp('---------------------------------------------');
disp(['There are ',num2str(Nfile), ' files in total']);

dataxAll = [];
dataRb2All = [];
dataRb2_bkg_All = [];
dataNcyc = [];

xdataAll = [];
ydataAll = [];
yerr_dataAll = [];

for i = 1:Nfile
    filename4load = [fnames(i).folder '/' fnames(i).name];
    load(filename4load);
    %     disp(['=> The ', num2str(i), 'th file Loaded.']);
    dataxi = xVarList1+foffset;                     %[GHz]
    dataRb2i = xVarCounts_Rb_2;             %ion counts of Rb2
    dataRb2i_bkg = xVarCountsBkgd_Rb_2;     %bkgnd ion counts of Rb2
    dataNcyci = cycleNumxVar;               %cycle number
    dataxyi = [dataxi dataRb2i dataRb2i_bkg dataNcyci];
    dataAll{i} = dataxyi;
    
    dataxAll = [dataxAll dataxi'];
    dataRb2All = [dataRb2All dataRb2i'];
    dataRb2_bkg_All = [dataRb2_bkg_All dataRb2i_bkg'];
    dataNcyc = [dataNcyc dataNcyci'];
end
disp('---------------------------------------------');


%%%------plot-------
f0 = 445000;

markersize = 6;
yplotmin = -0.5;
yplotmax = 30;
Xplot = dataxAll-f0;
Yplot = (dataRb2All-dataRb2_bkg_All)./dataNcyc;
Yplot_err = sqrt(dataRb2All+dataRb2_bkg_All)./dataNcyc;

% h2 = figure('units','normalized','WindowStyle','docked');

% % Yplot = (dataRb2All-dataRb2_bkg_All)./dataNcyc;
% % Yplot_err = sqrt(dataRb2All+dataRb2_bkg_All)./dataNcyc;

% errorbar(Xplot, Yplot, Yplot_err, 'ok');
% % ylim([yplotmin, yplotmax]);
% xlabel(['f - ', num2str(f0),' (GHz)']);
% ylabel('Ion count per cycle');
% title('All Rb_2 REMPI data with the rotation transitions of B^1\Pi_u(v''=4) - X^1\Sigma_g(v"=0,J")');
fcorrection = -0.17 - foffset;         %[GHz]
f_Q_branch = fcorrection+[...
    445000-1.867 445001.75 445005.19 445008.44 445011.52 445014.41 445017.12 ...
    445019.65 445022.01 445024.18-0.03 445026.17-0.03 445027.97-0.03 445029.6-0.03 ...
    445031.05-0.03 445032.31-0.03 445033.4-0.03  445034.3-0.03  445035.03-0.03 445035.57-0.03 ...
    445035.93-0.03 445037.24-0.03]; %[GHz]
col = 8;                        %plot column
NplotQ = length(f_Q_branch);
row = ceil(NplotQ/col);          %plot row

N_list = [];
A_list = [];
A_list_err = [];
dx_list = [];
dx_list_err = [];

h1 = figure('position', [0, 500, 1000, 600]);
subplot(2, 1, 1);
for j = 8:-1:2
    disp(' ');
    disp(['j = ', num2str(21-j)]);
    
    row = ceil(NplotQ/col);          %plot row
    fQc = f_Q_branch(j);
    dfmax = 0.3;                      %[GHz]
    dfmin = 0.3;
    if j==3
        dfmax = 0.3;                      %[GHz]
        dfmin = 0.18;
    end
    if j==5
        dfmax = 0.18;                      %[GHz]
        dfmin = 0.3;
    end
    xfitdata = Xplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata = Yplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata_err = Yplot_err(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    if 21-j==13
        xfitdata6 = xfitdata;
        yfitdata6 = yfitdata;
        yfitdata6_err = yfitdata_err;
    end
    
    %---fit data to a model----
    xc = fQc-f0;
    plot([xc,xc],[-1,40],'--b','LineWidth',0.5);
    hold on;
    
    Aguess = max(yfitdata);
    dxguess = 0.03;
    if mod(21-j,2) ==0
        fo = fitoptions('Method','NonlinearLeastSquares',...
            'Lower',[0, 0],...
            'Upper',[50, 0.07],...
            'StartPoint',[Aguess, dxguess]);
        ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','problem','xc','options',fo);
        [curvex,gofx] = fit(xfitdata', yfitdata', ft,'problem', xc);
    else
        fo = fitoptions('Method','NonlinearLeastSquares',...
            'Lower',[0, 0, xc-0.1],...
            'Upper',[50, 0.07, xc+0.1],...
            'StartPoint',[Aguess, xc, dxguess]);
        ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','options',fo);
        [curvex,gofx] = fit(xfitdata', yfitdata', ft)
    end
    N_list = [N_list 21-j];
    A_list = [A_list curvex.A];
    ci = confint(curvex);
    A_list_err = [A_list_err (ci(2,1)-ci(1,1))/4];
    dx_list = [dx_list curvex.dx];
    dx_list_err = [dx_list_err (ci(2,2)-ci(1,2))/4];
    
    curvex1 = xc-0.6:0.001:xc+0.6;
    plot(curvex1, curvex(curvex1), 'k', 'LineWidth', 1);
    if 21-j==13
        curvex6 = curvex1;
        curvey6 = curvex(curvex6);
    end
    legend off
    hold on
    
    if mod(21-j,2) ==1
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'MarkerFaceColor', [47 108 167]/(47+108+167)*1.9, 'CapSize', 0, 'LineWidth', 1);
    else
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'CapSize', 0, 'LineWidth', 1);
    end
    
    xdataAll = [xdataAll xfitdata];
    ydataAll = [ydataAll yfitdata];
    yerr_dataAll=[yerr_dataAll yfitdata_err];
    
    
    ax = gca;
    ax.FontSize = 13;
    xlabel(['f - ', num2str(f0),' (GHz)']);
    ylabel('Ion count per cycle');
    xlim([0, 21]);
    xticks(0:2:21);
    ylim([yplotmin, yplotmax]);
    
end

% h2 = figure('position', [100, 500, 1100, 700]);
subplot(2, 1, 2);
for j = 9:NplotQ
    disp(' ');
    disp(['j = ', num2str(21-j)]);
    row = ceil(NplotQ/col);          %plot row
    fQc = f_Q_branch(j);
    dfmax = 0.3;                      %[GHz]
    dfmin = 0.3;
    if j==10
        dfmax = 0.10;                      %[GHz]
        dfmin = 0.3;
    end
    if j==15
        dfmax = 0.13;                      %[GHz]
        dfmin = 0.3;
    end
    if j==19
        dfmax = 0.2;                      %[GHz]
        dfmin = 0.;
    end
    if j==20
        dfmax = 0.3;                      %[GHz]
        dfmin = 0.18;
    end
    xfitdata = Xplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata = Yplot(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    yfitdata_err = Yplot_err(Xplot<=fQc-f0+dfmax&Xplot>=fQc-f0-dfmin);
    
    
    if 21-j==2
        yfitdata = yfitdata-19.52*exp(-(xfitdata-35.30-0.08).^2/2/0.0408^2);
    end
    if 21-j==6
        yfitdata = yfitdata-16.41*exp(-(xfitdata-32.24-0.08).^2/2/0.03965^2);
    end
    
    %---fit data to a model----
    xc = fQc-f0;
    plot([xc,xc],[-1,40],'--b','LineWidth',0.5);
    hold on;
    Aguess = max(yfitdata);
    dxguess = 0.03;
    if mod(21-j,2) ==0
        fo = fitoptions('Method','NonlinearLeastSquares',...
            'Lower',[0, 0.04],...
            'Upper',[50, 0.07],...
            'StartPoint',[Aguess, dxguess]);
        ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','problem','xc','options',fo);
        [curvex,gofx] = fit(xfitdata', yfitdata', ft,'problem', xc);
    else
        fo = fitoptions('Method','NonlinearLeastSquares',...
            'Lower',[0, 0.04, xc-0.1],...
            'Upper',[50, 0.07, xc+0.1],...
            'StartPoint',[Aguess, xc, dxguess]);
        ft = fittype('A.*exp(-(x-xc).^2/2/dx^2)','options',fo);
        [curvex,gofx] = fit(xfitdata', yfitdata', ft)
    end
    N_list = [N_list 21-j];
    A_list = [A_list curvex.A];
    ci = confint(curvex);
    A_list_err = [A_list_err (ci(2,1)-ci(1,1))/4];
    dx_list = [dx_list curvex.dx];
    dx_list_err = [dx_list_err (ci(2,2)-ci(1,2))/4];
    
    
    curvex1 = xc-0.3:0.001:xc+0.3;
    plot(curvex1, curvex(curvex1), 'k', 'LineWidth', 1);
    
    hold on
    legend off
    if mod(21-j,2) ==1
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'MarkerFaceColor', [47 108 167]/(47+108+167)*1.9, 'CapSize', 0, 'LineWidth', 1);
    else
        errorbar(xfitdata, yfitdata, yfitdata_err, 'ok', 'MarkerSize', markersize,...
            'CapSize', 0, 'LineWidth', 1);
    end
    
    xdataAll = [xdataAll xfitdata];
    ydataAll = [ydataAll yfitdata];
    yerr_dataAll=[yerr_dataAll yfitdata_err];
    
    ax = gca;
    ax.FontSize = 13;
    xlabel(['f - ', num2str(f0),' (GHz)']);
    ylabel('Ion count per cycle');
    xlim([21, 40]);
    xticks(21:2:40);
    ylim([yplotmin, yplotmax]);
end


if 0
    axes('Position',[0.73 0.23 0.15 0.2]);
    box on;
    plot(curvex6, curvey6, 'k', 'LineWidth', 1);
    hold on
    errorbar(xfitdata6, yfitdata6, yfitdata6_err, 'ok', 'MarkerSize', markersize+2,...
        'MarkerFaceColor', [47 108 167]/(47+108+167)*1.9, 'CapSize', 0, 'LineWidth', 1);
    hold off
    xlim([19.30, 19.80]);
    ylim([yplotmin, 23]);
    ax = gca;
    ax.FontSize = 11;
end

h2 = figure('position', [0, 500, 700, 250]);
odd = mod(N_list,2);
even = ~odd;
bar_x = N_list;
bar_y = (A_list-0.26)*sqrt(2*pi).*dx_list/0.03;
bar_y_err = bar_y.*min(sqrt((A_list_err./A_list).^2),1);%+(dx_list_err./dx_list).^2
% bar(bar_x((bar_x<20)),bar_y((bar_x<20)), 'FaceColor', [0 1 1]*0.5, 'EdgeColor', [0 1 1]);
bar(bar_x((bar_x<21)),bar_y((bar_x<21)), 'FaceColor', [47 108 167]/(47+108+167)*1.3, 'EdgeColor', [47 108 167]/(47+108+167)*1.);
yticks([0:95/5:95]);
yticklabels({'0','0.2','0.4', '0.6','0.8','1'})
hold on
bar_y1=bar_y;
bar_y1((bar_x<20)&odd)=0;
bar(bar_x(bar_x<20),bar_y1(bar_x<20), 'FaceColor', [102 8 32]/(102+8+32)*0.7, 'EdgeColor', [102 8 32]/(102+8+32)*0.5);
% bar(bar_x(bar_x<20),bar_y1(bar_x<20), 'FaceColor', 'red', 'EdgeColor', [1 0 1]);
xticks([0:20]);
xlim([-1 20.5]);
ylim([0 95]);
errorbar(bar_x(bar_x<20),bar_y(bar_x<21),bar_y_err(bar_x<21), '.', 'MarkerSize', 0.5,...
    'CapSize', 3, 'LineWidth', 1, 'MarkerEdgeColor',[1 1 1]*0,...
    'Color',[1 1 1]*0);
hold off
% xlabel('Rotational state N');
% ylabel('Integrated ion count per cycle');
% title('K_2 REMPI');
figure(h2)

set(h1,'Units','Inches');
pos = get(h1,'Position');
set(h1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h1,'plots\K2REMPI_spectrum_30G_inset','-dpdf','-r0')

set(h2,'Units','Inches');
pos = get(h2,'Position');
set(h2,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h2,'plots\Rb2REMPI_bar_30G','-dpdf','-r0')